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TITLE QF THE INVENTION 

INTERPOLATION METHOD , APPARATUS FOR CARRYING OUT THE 
METHOD, AND CONTROL PROGRAM FOR IMPLEMENTING THE METHOD 

BACKGROUND QF THE INVENTION 

Field of the Invention 

The present invention relates to an interpolation 
method that is used in computational fluid dynamics for 
realizing the motion of a fluid on a computer, an 
apparatus for carrying out the method, and a program 
for implementing the method. 

Description of the Related Art 

In recent years, computational fluid dynamics that 
realizes the motion of a fluid on a computer has made 
remarkable progress. However, it is difficult to 
describe the motion of an object in a fluid, or the 
motion of two fluids or a multi-layer fluid, and 
therefore, various approaches have been made to this 
difficult theme. 

To compute a fluid flow around a stationary object, 
a grid can be generated according to the shape of the 
object, and it is a conventionally known method of the 
computation to utilize an unstructured grid or the like 
based on a finite element method or the like. When the 
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object is moved, a method of changing the grid 
according to the motion of the object can be employed. 
However, if the change of the grid involves large 
deformation, it is generally difficult to perform the 
5 computation. 

On the other hand, a difference equation that 
approximately realizes a fluid equation describing a 
fluid flow on a structured grid, by using duality 
between a grid (primary grid) and a dual grid known as 

10 a staggered grid, is highly reliable due to its 

stability. This leads to an idea of describing a 
volume of a substance in each cell as a function on a 
grid, and moving the function. One example of the 
method based on this idea utilizes the Volume-of -Fluid 

15 method (hereinafter referred to as "the VOF method") 
for computing a flow of an incompressible fluid. In 
this case, it is a key technique to move a function 
describing the volume of the substance of the fluid. 

Further, in the case of a compressible fluid, the 

2 0 computation of a flow thereof involves key problems of 
defining motions of various objects, and in particular, 
propagation of a shock wave and the like are very 
closely related with the above-mentioned theme of 
describing the motion of a substance or an interface. 

25 An initial-value problem concerning an equation 

(referred to as "an advection equation") describing 
such a motion is called the Riemann problem after a 
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paper by Riemann on a wave equation concerning a wave 
in the air "Ueber die Fortpf lanzung ebener Luftwellen 
von endlicher Schwingungsweite (Albhandlungen der 
Koniglichen Gesellschaft der Wissenschaf ten zu 
5 Gottingen, 8 (I860), 43-65". 

As numerical solutions to the Riemann problem, 
there have been proposed various concepts and schemes, 
independent of each other or dependent on each other, 
including a high-order accuracy upwind scheme called 

10 the Total Variation Diminishing scheme (hereinafter 

referred to as "the TVD scheme"), the Essentially Non- 
Oscillatory scheme (hereinafter referred to as "the ENO 
scheme") which introduces a limiter function, the 
Monotone Upstream Centered Scheme for Conservation Laws 

15 (hereinafter referred to as "the MUSCL scheme", and the 
Positivity-Preserving scheme (hereinafter referred to 
as "the PP scheme". 

One of equations dealt with in the present 
invention is a one-dimensional partial differential 

2 0 equation defined by: 

dF(x,t)/<? t =dC(x, t)F(x, t) / dx ...(1) 
To perform computation based on the equation, we 
will now consider the problem of converting this 
equation into a difference equation. In the following, 

25 arguments X, T of a function F(X,T) represent a 

discrete grid point in space and a discrete time point, 
respectively . 
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To solve the equation (1), there have been 
proposed various computing methods in the computational 
fluid dynamics as described above. The TVD scheme is 
to perform the computation under a TVD condition that a 
5 TVD value TV(F[-,n]) defined by: 

TV(F[-,n] )=Sj |F[j+l,n]-F[j,n] | .-.(2) 
does not increase with time evolution, the TVD 
condition being defined by the following inequality: 

TV(F[ • ,n+l] )<TV(F[- ,n] ) ...(3) 
10 The MUSCL scheme, in a broader sense, is a 

computational method composed of data reconstruction 
and an approximate Riemann solution. 

The ENO scheme is an upwind method with a limiter 
function in which the MUSCL scheme and the TVD scheme 
15 are combined to eliminate vibration occurring when 
there are discontinuities in the above-described 
function. 

Further, for data reconstruction, there are 
conventionally employed a method of interpolating the 

20 above-described function with a piecewise constant 

function (which gives first-order accuracy) , a method 
of interpolating the above-described function with a 
piecewise linear function (which gives second-order 
accuracy) , and a method of interpolating the above- 

25 described function by a piecewise parabolic function 
(which gives third-order accuracy) . 

The present invention relates to improvement in 
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the second-order accuracy of interpolation. 

Now, a description will be given of a conventional 
second-order ENO scheme. 

First, a description will be given of a method of 
5 numerically solving the above equation (1) by finding a 
solution thereto by finite difference approximation, 
assuming that there is a piecewise linear function: 

H[I,J] (x)=F(I, J)+G(I, J)x ...(4) 

If the function C( 1+1/2, J) is positive, there is 
10 given the following definition: 

Q (1 + 1/2, J) :=C( 1+1/2, J)H[I, J] ( Ax/2-C (I, J) At/2) 

...(5) 

If the function is not positive, there is given 
the following definition: 
15 Q(I+l/2, J) 

:=C(I+l/2, J)H[I + 1, J] (-Ax/2-C (1 + 1, J) At/2) ...(6) 
The above equation (1) can be approximately solved 
by obtaining, from the above definitions, the following 
equation: 
20 F(I,J+1) 

=F (I, J) + A t (Q (1+1/2 , J) -Q (1-1/2, J) ) / Ax ...(7) 
On the other hand, a slope function G(I,J) 
representative of the slope of the piecewise linear 
function of the equation (4) is defined by the 
25 following equation: 

G(I, J)=AF(I, J) / Ax ...(8) 

Now, there are various methods of defining 



AF(I,J), as described hereinafter, and typical three 
of them will be described here. 

First, there are introduced the following 
definitions : 

Sp(I,J) : =F ( 1+1 , J) -F ( I , J) 

Sm(I,J) :=F(I,J)-F(I-1,J) ...(9) 

A TVD limiter function in one dimension called the 
minmod function gives AF(I,J) defined by: 
A F ( I , J ) = ( sgn ( Sp ( I , J ) ) +sgn (Sm ( I , J) ) ) 
Min( |Sp(I,J) | , |Sm(I, J) | ) /2 ...(10) 

A TVD limiter function in one dimension called the 
averaging limiter function gives Af(I,J) defined by: 
A F ( I , J ) = ( sgn ( Sp ( I , J ) ) +sgn (Sm (I , J) ) ) 
Min( | Sp (I f J") +Sm(I, J) | /2 , 2 | Sp (I, J) | , 2 \ Sm(I, J) | ) /2 

...(11) 

A TVD limiter function in one dimension called the 
superbee limiter function gives AF(I,J) defined by: 

A F ( I , J ) = ( sgn ( Sp ( I , J ) ) +sgn(Sm(I / J) ) ) 

Min(Max( | Sp (I , J) | , | Sm ( I , J) | ) ) , 

2 |Sp(I, J) | ,2 |Sm(I, J) | ) /2 ...(12) 

In the above equations (10) to (12), sgn( ) 
represents a sign function. That is, when the argument 
thereof is positive i.e. of a plus sign, the function 
gives a value of 1, while the same is negative i.e. of 
a minus sign, the function gives a value of -1. 

These methods are described e.g. in "A. Suresh, 
"Positivity-Preserving Schemes in Multidimensions" SIAM 



J. Sci. Comput. 22, 1184-1198(2000)". 

Using these interpolation methods, a function 
describing a square shape was subjected to advection 
calculation using the above equation (7) . FIGS. 7 to 9 
5 show results of the calculation. More specifically, in 
FIGS. 7 to 9, shapes are depicted which are given by 
the function after 100 steps under periodic boundary 
conditions, using the minmod limiter function, the 
averaging limiter function, and the superbee limiter 
10 function, together with the initial shape for 

comparison. As is clear from FIGS. 7 to 9, the square 
shape as the initial shape is drastically lost into 
• dull shapes. 

Next, a description will be given of an algorithm 
15 for a higher dimension obtained by extending the 

algorithm for one dimension described heretofore. A 
method employed in the algorithm for a higher dimension 
described here is disclosed in the above-mentioned 
literature. 

20 Suppose that F(I,J,T) represents a function of a 

point on a two-dimensional grid at a given time T, 
wherein I and J are integers and cooperatively 
represent a grid point on the two-dimensional grid. 

For the function, there are defined the following 

25 values: 

Vmin=Min<F(I-l, J+1,T) -F(I, J,T) , F(I, J+1,T) -F(I, J,T) , 
F(I+1, J+1,T)-F(I, J,T) ,F(I-1, J,T)-F(I, J,T) , - 8 , 
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F(I+1, J,T)-F(I, J,T) ,F(I-1, J-1,T)-F(I, J,T) , 
F(I, J-1,T)-F(I, J,T) ,F(I+1, J-1,T)-F(I, J,T) ) 
Vmax=Max(F(I-l, J+1,T) -F(I, J,T) ,F(I, J+1,T) -F(I, J,T) , 
F(I+1, J+1,T) -F(I, J,T) ,F(I-1, J,T)-F(I, J,T) , £ , 
F(I+1,J,T)-F(I,J,T) ,F(M,J-1,T)-F(I,J,T) , 
F(I, J-l.T) -F(I, J,T) ,F(I+1,J-1,T) -F(I, J,T) ) ...(13) 
Further, there are defined: 
Sx(I, J,T)=(F(I+1, J,T) -F(I-1, J,T) ) /2 
Sy(I, J,T)=(F(I, J+1,T)-F(I, J-1,T) )/2 ...(14) 
V=2 (Min( |Vmin| , |Vmax| ) ) / ( |Sx(I, J,T) | 
+ |Sx(I, J,T) | + £ ) ...(15) 
At this time, after defining: 
Gx(I, J,T)=Min(l,V) Sx(I, J,T) / Ax 
Gy(I, J,T) =Min(l,V) Sy ( I , J, T) / Ay ...(16) 
an interpolation function H(I, J,T) (x,y) for a grid 
point (l,j) is defined as: 

H(I,J,T) (x,y)=F(I, J,T)+Gx(I, J,T)x+Gy(I, J,T)y ...(17) 
By using the method, with respect to initial 
values of a function shown in FIG. 10A, advection 
calculation of the function can be carried, a result of 
which is shown in FIG. 10B. 

However, the conventional interpolation method 
suffers from the following problem: 

In the Volume-of -Fluid (VOF) method applied to a 
two-phase incompressible fluid, how to cause time 
evolution of a shape-describing function descriptive of 
a shape as a result of the advection calculation while 
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preserving the shape described by the function is a 
matter of great interest, but the conventional 
interpolation method tends to cause progressive loss of 
the initial shape into a dull shape, with execution of 
5 computational steps. 

In view of this problem, a computational algorithm 
enabling the shape to be preserved more effectively is 
indispensable, and to this end, for an algorithm for 
one dimension, it is necessary to replace the 
10 interpolation function of the equation (4) by an 
improved one, and for an algorithm for a higher 
dimension, it is necessary to replace the interpolation 
function of the equation (17) by an improved one. 

15 SUMMA R Y O F THE INVENTION 

The present invention has been made in view of the 
above problems of the prior art, and it is an object of 
the present invention is to provide an interpolation 

2 0 method for the Volume-of -Fluid (VOF) method applied to 
a two-phase incompressible fluid, in particular, which 
makes it possible to cause time evolution of a shape- 
describing function based on the Volume-of -Fluid (VOF) 
method while preserving a sharpness of a shape 

2 5 described by the function, an apparatus for carrying 

out the method, and a control program for implementing 
the method. 
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To attain the above object, in a first aspect of 
the present invention, there is provided an 
interpolation method of defining a function F on a one- 
dimensional structured grid formed on a one-dimensional 
5 real region, the function being defined through 

definition of a value thereof at a center of each cell 
within the one-dimensional structured grid, as an 
interpolation function H, the method comprising the 
steps of setting, with respect to a cell of interest on 

10 the one-dimensional structured grid, a slope to zero if 
a forward difference and a backward difference of the 
function f have different signs, and to a value twice 
as large as a smaller one of absolute values of the 
forward difference and the backward difference if the 

15 forward difference and the backward difference have the 
same sign, and defining the function F on a partial 
region of the one-dimensional real region determined by 
the cell of interest, by a linear function having a 
value of F0 at a center of the cell of interest and the 

20 slope. 

With the arrangement of the interpolation method 
according to the first aspect of the present invention, 
it is possible to realize an interpolation algorithm 
for the Volume-of -Fluid (VOF) method applied to a two- 
25 phase incompressible fluid, in particular, which makes 
it possible to cause time evolution of a shape- 
describing function based on the Volume-of -Fluid (VOF) 
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method while preserving a sharpness of a shape 
described by the function, and the function gives very- 
excellent results of describing the shape. 

Preferably, the interpolation method is applied to 
5 as a numerical solution of an advection-type 
differential equation . 

To attain the above object, in a second aspect of 
the present invention, there is provided an 
interpolation method of defining a function F on a two- 

10 dimensional structured grid formed on a two-dimensional 
real region, the function being defined through 
definition of a value thereof at a center of each cell 
within the two-dimensional structured grid, as an 
interpolation function H, the method comprising the 

15 steps of setting a cell of interest to a cell A, the 

cell A having a first side extending in an x direction, 
a second side extending in the x direction and being 
opposite to the first side, a third side extending in a 
y direction, and a fourth side extending in the y 

20 direction and being opposite to the third side, 

defining values twice as large as values of one-sided 
differences of the function F between a center of the 
cell A and respective centers of four cells adjacent to 
the cell A on the first side, the second side, the 

25 third side, and the fourth side, as a first-sided 

difference (DFxmin) , a second-sided difference (DFxmax) , 
a third-sided difference (DFymin) , and a fourth-sided 
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difference (DFymax) , respectively, and setting an x- 
direction difference to zero if the first-sided 
difference and the second-sided difference have 
different signs, and to a smaller one of absolute 
5 values of the first-sided difference and the second- 
sided difference if the first-sided difference and the 
second-sided difference have the same sign, and a in- 
direction difference to zero if the third-sided 
difference and the fourth-sided difference have 

10 different signs, and to a smaller one of absolute 

values of the third-sided difference and the fourth- 
sided difference if the third-sided differences and the 
fourth-sided difference have the same sign, forming an 
interpolation function candidate which has slopes 

15 defined by the x-direction difference and the y- 

direction difference, respectively, and has a value F0 
of the function F at the center of the cell A and 
setting this candidate to a plane B, modifying, if a 
value of the plane B at a first grid point at a 

20 location of intersection of the first side and the 

third side of the cell A is larger than a value of the 
center of the cell A, the plane B by multiplying the x- 
direction difference and the y-direction difference by 
a largest constant not more than 1 such that the value 

25 of the plane B at the first grid point does not exceed 
any of values of the function F at respective centers 
of three cells having the first grid point in common 
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except for the cell A, modifying, if the value of the 
plane B at the first grid point is smaller than the 
value of the center of the cell A, the plane B by 
multiplying the x-direction difference and the y- 
5 direction difference by a largest constant not more 

than 1 such that the value of the plane B at the first 
grid point does not fall below any of the values of the 
function F at the respective centers of the three cells 
having the grid point in common except for the cell A, 

10 and carrying out, on the plane B thus obtained, the 
same operation as carried out as to the first grid 
point, as to a second grid point at a location of 
intersection of the first side and the fourth side of 
the cell A, a third grid point at a location of 

15 intersection of the second side and the third side, and 
a fourth grid point at a location of intersection of 
the second side and the fourth side, to thereby change 
the slope of the plane B, and define the resulting 
plane as the interpolation function within the cell A. 

2 0 With the arrangement of the interpolation method 

according to the second aspect of the present invention, 
it is possible to obtain the same advantageous effects 
as provided by the first aspect of the present 
invention. 

25 To attain the above object, in a third aspect of 

the present invention, there is provided an 
interpolation method of defining a function F on a 



three-dimensional structured grid formed on a three- 
dimensional real region, the function being defined 
through definition of a value thereof at a center of 
each cell within the three-dimensional structured grid, 
5 as an interpolation function H, the method comprising 
the steps of setting a cell of interest to a cell A, 
the cell A having a first side extending in an x 
direction, a second side extending in the x direction 
and being opposite to the first side, a third side 

10 extending in a y direction, a fourth side extending in 
the y direction and being opposite to the third side, a 
fifth side in a 2 direction, and a sixth side extending 
in the z direction and being opposite to the fifth side, 
defining values twice as large as values of one-sided 

15 differences of the function F between a center of the 
cell A and respective centers of six cells adjacent to 
the cell A on the first side, the second side, the 
third side, the fourth side, the fifth side, and the 
sixth side, as a first-sided difference (DFxmin) , a 

20 second-sided difference (DFxmax) , a third-sided 

difference (DFymin) , a fourth-sided difference (DFymax) , 
a fifth-sided difference (DFzmax) , and a sixth-sided 
difference (DFzmin) , respectively, and setting an x- 
direction difference to zero if the first-sided 

25 difference and the second-sided difference have 

different signs, and to a smaller one of absolute 
values of the first-sided difference and the second- 
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grid point does not exceed any of the values of the 
function F at respective centers of seven cells having 
the first grid point in common except for the cell A, 
modifying, if the value of the plane B at the first 
5 grid point is smaller than the value of the center of 
the cell A, the plane B by multiplying the x-direction 
difference, the y-direction difference, and the z- 
direction difference, by a largest constant not more 
than 1 such that the value of the plane B at the first 

10 grid point does hot fall below any of the values of the 
function F at the respective centers of the seven cells 
having the grid point in common except for the cell A, 
and carrying out, on the plane B thus obtained, the 
same operation as carried out as to the first grid 

15 point, on a second grid point at a location of 

intersection of the first side, the third side, and the 
sixth side, a third grid point at a location of 
intersection of the first side, the fourth side, and 
the fifth side, a fourth grid point at a location of 

2 0 intersection of the first side, the fourth side, and 
the sixth side, a fifth grid point at a location of 
intersection of the second side, the third side, and 
the fifth side, and a sixth grid point at a location of 
intersection of the second side, the third side, and 

25 the sixth side, a seventh grid point at a location of 
intersection of the second side, the fourth side, and 
the fifth side, and an eighth grid point at a location 



17 

of intersection of the second side, the fourth side, 
and the sixth side, to thereby change the slope of the 
plane, and define the resulting plane as the 
interpolation function within the cell A. 
5 With the arrangement of the interpolation method 

according to the. third aspect of the present invention, 
it is possible to obtain the same advantageous effects 
as provided by the first aspect of the present 
invention. 

10 To attain the above object, in a fourth aspect of 

the present invention, there is provided an apparatus 
for carrying out an interpolation method of defining a 
function F on a one -dimensional structured grid formed 
on a one-dimensional real region, the function being 

15 defined through definition of a value thereof at a 
center of each cell within the one-dimensional 
structured grid, as an interpolation function H, the 
apparatus comprising a setting device that sets, with 
respect to a cell of interest on the one-dimensional 

20 structured grid, a slope to zero if a forward 

difference and a backward difference of the function F 
have different signs, and to a value twice as large as 
a smaller one of absolute values of the forward 
difference and the backward difference, and a 

25 definition device that defines the function F on a 
partial region of the one -dimensional real region 
determined by the cell of interest, by a linear 
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function having a value of FO at a center of the cell 
of interest and the slope. 

With the arrangement of the apparatus according to 
the fourth aspect of the present invention, it is 
5 possible to obtain the same advantageous effects as 
provided by the first aspect of the present invention. 

To attain the above object, in a fifth aspect of 
the present invention, there is provided an apparatus 
for carrying out an interpolation method of defining a 

10 function F on a two-dimensional structured grid formed 
on a two-dimensional real region, the function being 
defined through definition of a value thereof at a 
center of each cell within the two-dimensional 
structured grid, as an interpolation function H, the 

15 apparatus comprising a cell-setting device that sets a 
cell of interest to a cell A, the cell A having a first 
side extending in an x direction, a second side 
extending in the x direction and being opposite to the 
first side, a third side extending in a y direction, 

20 and a fourth side extending in the y direction and 

being opposite to the third side, a difference-setting 
device that defines values twice as large as values of 
one-sided differences of the function F between a 
center of the cell A and respective centers of four 

25 cells adjacent to the cell A on the first side, the 

second side, the third side, and the fourth side, as a 
first-sided difference (DFxmin) , a second-sided 



19 

difference (DFxmax) , a third-sided difference (DFymin) , 
and a fourth-sided difference (DFymax) , respectively, 
and sets an x-direction difference to zero if the 
first-sided difference and the second-sided difference 
5 have different signs, and to a smaller one of absolute 
values of the first-sided difference and the second- 
sided difference if the first-sided difference and the 
second-sided difference have the same sign, and a y- 
direction difference to zero if the third-sided 

10 difference and the fourth-sided difference have 

different signs, and to a smaller one of absolute 
values of the third-sided difference and the fourth- 
sided difference if the third-sided differences and the 
fourth-sided difference have the same sign, an 

15 interpolation function candidate -forming device that 
forms an interpolation function candidate which has 
slopes defined by the x-direction difference and the y- 
direction difference, respectively, and has a value F0 
of the function F at the center of the cell A and 

20 setting this candidate to a plane B, a first 

modification device that modifies, if a value of the 
plane B at a first grid point at a location of 
intersection of the first side and the third side of 
the cell A is larger than a value of the center of the 

25 cell A, the plane B by multiplying the x-direction 

difference and the y-direction difference by a largest 
constant not more than 1 such that the value of the 
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plane B at the first grid point does not exceed any of 
values of the function F at respective centers of three 
cells having the first grid point in common except for 
the cell A, a second modification device that modifies, 
5 if the value of the plane B at the first grid point is 
smaller than the value of the center of the cell A, the 
plane B by multiplying the x-direction difference and 
the y-direction difference by a largest constant not 
more than 1 such that the value of the plane B at the 

10 first grid point does not fall below any of the values 
of the function F at the respective centers of the 
three cells having the grid point in common except for 
the cell A, and a definition device that carries out, 
on the plane B thus obtained, the same operation as 

15 carried out as to the first grid point, as to a second 
grid point at a location of intersection of the first 
side and the fourth side of the cell A, a third grid 
point at a location of intersection of the second side 
and the third side, and a fourth grid point at a 

20 location of intersection of the second side and the 

fourth side, to thereby change the slope of the plane B, 
and define the resulting plane as the interpolation 
function within the cell A. 

With the arrangement of the apparatus according to 

25 the fifth aspect of the present invention, it is 

possible to obtain the same advantageous effects as 
provided by the first aspect of the present invention. 
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To attain the above object, in a sixth aspect of 
the present invention, there is provided an apparatus 
for carrying out an interpolation method of defining a 
function F on a three-dimensional structured grid 
5 formed on a three-dimensional real region, the function 
being defined through definition of a value thereof at 
a center of each cell within the three-dimensional 
structured grid, as an interpolation function H, the 
apparatus comprising a cell-setting device that sets a 

10 cell of interest to a cell A, the cell A having a first 
side extending in an x direction, a second side 
extending in the x direction and being opposite to the 
first side, a third side extending in a y direction, a 
fourth side extending in the y direction and being 

15 opposite to the third side, a fifth side in a z 
direction, and a sixth side extending in the z 
direction and being opposite to the fifth side, a 
difference-setting device that defines values twice as 
large as values of one-sided differences of the 

20 function F between a center of the cell A and 

respective centers of six cells adjacent to the cell A 
on the first side, the second side, the third side, the 
fourth side, the fifth side, and the sixth side, as a 
first-sided difference (DFxmin) , a second-sided 

25 difference (DFxmax) , a third- sided difference (DFymin) , 
a fourth-sided difference (DFymax) , a fifth-sided 
difference (DFzmax) , and a sixth-sided difference 
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(DFzmin) , respectively, and sets an x-direction 
difference to zero if the first-sided difference and 
the second-sided difference have different signs, and 
to a smaller one of absolute values of the first-sided 
5 difference and the second-sided difference if the 

first-sided difference and the second-sided difference 
have the same sign, a y-direction difference to zero if 
the third-sided difference and the fourth-sided 
difference have different signs, and to a smaller one 

10 of absolute values of the third-sided difference and 
the fourth-sided difference if the third-sided 
differences and the fourth-sided difference have the 
same sign, and a z-direction difference to zero if the 
fifth-sided difference and the sixth-sided difference 

15 have different signs, and to a smaller one of absolute 
values of the fifth-sided difference and the sixth- 
sided difference if the fifth-sided difference and the 
sixth-sided difference have the same sign, an 
interpolation function candidate- forming device that 

20 forms an interpolation function candidate which has 
slopes defined by the x-direction difference, the y- 
direction difference, and the z-direction difference, 
respectively, and has a value F0 of the function F at 
the center of the cell A, and setting this candidate to 

25 a plane B, a first modification device that modifies, 
if a value of the plane B at a first grid point at a 
location of intersection of the first side, the third 
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side, and the fifth side of the cell A is larger than 
the value of the center of the cell A, the plane B by 
multiplying the x-direction difference, the y-direction 
difference, and the z-direction difference, by a 
5 largest constant not more than 1 such that the value of 
the plane B at the first grid point does not exceed any 
of the values of the function F at respective centers 
of seven cells having the first grid point in common 
except for the cell A, a second modification device 

10 that modifies, if the value of the plane B at the first 
grid point is smaller than the value of the center of 
the cell A, the plane B by multiplying the x-direction 
difference, the y-direction difference, and the z- 
direction difference, by a largest constant not more 

15 than 1 such that the value of the plane B at the first 
grid point does not fall below any of the values of the 
function F at the respective centers of the seven cells 
having the grid point in common except for the cell A, 
and a definition device that carries out, on the plane 

20 B thus obtained, the same operation as carried out as 
to the first grid point, on a second grid point at a 
location of intersection of the first side, the third 
side, and the sixth side, a third grid point at a 
location of intersection of the first side, the fourth 

25 side, and the fifth side, a fourth grid point at a 

location of intersection of the first side, the fourth 
side, and the sixth side, a fifth grid point at a 
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location of intersection of the second side, the third 
side, and the fifth side, and a sixth grid point at a 
location of intersection of the second side, the third 
side, and the sixth side, a seventh grid point at a 
5 location of intersection of the second side, the fourth 
side, and the fifth side, and an eighth grid point at a 
location of intersection of the second side, the fourth 
side, and the sixth side, to thereby change the slope 
of the plane , and define the resulting plane as the 

10 interpolation function within the cell A. 

With the arrangement of the apparatus according to 
the sixth aspect of the present invention, it is 
possible to obtain the same advantageous effects as 
provided by the first aspect of the present invention. 

15 To attain the above object, in a seventh aspect of 

the present invention, there is provided a control 
program for causing a computer to execute the 
interpolation method as recited hereinabove as the 
interpolation method according to the first aspect of 

20 the present invention. 

To attain the above object, in an eighth aspect of 
the present invention, there is provided a control 
program for causing a computer to execute the 
interpolation method as recited hereinabove as the 

25 interpolation method according to the second aspect of 
the present invention. 

To attain the above object, in a ninth aspect of 
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the present invention, there is provided a control 
program for causing a computer to execute the 
interpolation method as recited hereinabove as the 
interpolation method according to the third aspect of 
5 the present invention. 

With the arrangement of each of the control 
programs according to the seventh to ninth aspects of 
the present invention, it is possible to obtain the 
same advantageous effects as provided by the first 
10 aspect of the present invention. 

The above and other objects, features, and 
advantages of the invention will become more apparent 
from the following detailed description taken in 
conjunction with the accompanying drawings. 

15 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIG . 1 is a flowchart showing a process of 
advection calculation using an interpolation method 
2 0 according to a first embodiment of the present 
invention; 

FIG. 2 is a diagram showing results of the 
advection calculation using the interpolation method 
illustrated in FIG. 1; 
25 FIG. 3 is a perspective view of a computer system 

for executing the advection calculation using the 
interpolation method illustrated in FIG. 1; 
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FIG. 4 is a block diagram showing the arrangement 
of essential components of the computer system using 
the interpolation method illustrated in FIG. 1; 

FIG. 5 is a flowchart showing a process of 
5 advection calculation using an interpolation method 
according to a second embodiment of the present 
invention; 

FIGS. 6A and 6B are diagrams showing results of 
the advection calculation using the interpolation 
10 method according to the second embodiment and an 

interpolation method according to a third embodiment of 
the present invention; 

FIG. 7 is a diagram showing results of advection 
calculation using a conventional minmod limiter 
15 function; 

FIG. 8 is a diagram showing results of advection 
calculation using a conventional averaging limiter 
function; 

FIG. 9 is a diagram showing results of advection 
20 calculation using a conventional superbee limiter 
function; and 

FIGS. 10A and 10b are diagrams showing results of 
advection calculation using a conventional 
interpolation method (for a higher dimension) . 

25 



DETAILED DESCRIPTION QF THE PREFERRED EMBODIMENTS 
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The present invention will now be described in 
detail below with reference to the accompanying 
drawings showing preferred embodiments thereof . In the 
drawings, elements and parts which are identical 
5 throughout the views are designated by identical 

reference numeral, and duplicate description thereof is 
omitted. 

First, an interpolation method according to a 
first embodiment of the invention will be described. 

10 The method according to the present embodiment 

interpolates a function by a piecewise linear function, 
for data reconstruction. That is, as a numerical 
solution to an advection equation, the present 
interpolation method corresponds to a solution of the 

15 second-order accuracy. 

In the present embodiment, as a problem for one 
dimension, we will consider the problem of converting a 
one-dimensional partial differential equation defined 
by: 

20 dF(x,t)/d t = dC(x,t)F(x,t)/dx ...(18) 

into a difference equation to perform computation based 
on the equation. In the following, arguments X, T of a 
function F(X,T) represent a discrete grid point in 
space and a discrete time point, respectively. 

25 First, it is assumed that the interpolation 

function H [1/ J] (x) according to the present embodiment 
has been already obtained as a piecewise linear 
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function as follows: 

H[I, J] (x)=F(I / J)+G(I, J)x ...(19) 

With this assumption, a description will be given 
of a method of numerically solving the above equation 
5 (18) by finding a solution thereto by finite difference 
approximation, A method of determining the coefficient 
G(I,J) of the interpolation function H[I,J](x) of the 
present embodiment will be described hereinafter. As 
described later, the method of determining the 

10 coefficient G(I,J) characterizes the interpolation 

algorithm according to the present embodiment. It is 
essential to the present embodiment that the 
coefficient G(I,J) is determined by using neighboring 
information of geometric properties of the grid. 

15 In short, when there are provided a one- 

dimensional structured grid formed on a one-dimensional 
real region, and a function F(X,T) on the one- 
dimensional structured grid, the function being defined 
through definition of a value thereof at the center of 

20 each cell within the one-dimensional structured grid, 
the function F(X,T) is defined as an interpolation 
function H[I,J](x) on the one-dimensional real region 
by using the function F(X,T) and the geometry of the 
grid. In other words, while the function F defines 

25 mapping by discrete points alone, the function H 

extends the definition range of mapping over the entire 
X axis . 
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If the function C (1+1/2 , J) is positive, there is 
given the following definition: 

Q (1 + 1/2 , J) :=C(I+1/2,J)H[I,J] (Ax/2-C(I,J) At/2) 

..(20) 

If the same is not positive, there is given the 
following definition: 

Q ( 1+1/2 , J) : =C (I+1/2,J)H[I + 1,J] ( - Ax/2-C (I , J) At/2) 

...(21) 

From the above definitions, there is obtained the 
following equation: 

F(I,J+l)=F(I,J)+At(Q(I+l/2,J) -Q (1-1/2 , J) ) / Ax 

...(22) 

On the other hand, the slope function G(I,J) 
representative of the slope of the piecewise linear 
function is defined in the following manner: 

First, there are introduced the following 
definitions : 

Sp ( I , J) :=F(I+1, J)-F(I, J) 

Sm(I,J) : =F ( I , J) -F ( 1-1 , J) ...(23) 

Further, there is introduced a limiter function in 
one dimension for the present embodiment, defined by: 

A F ( I , J ) = ( sgn ( Sp { I , J ) ) +sgn ( Sm ( I , J) ) ) 

Min ( | Sp ( I , J) | , |Sm(I, J) | ) ...(24) 

In this equation, sgn( ) represents a sign 
function. That is, when the argument thereof is 
positive or of a plus sign, the function gives a value 
of 1, while the same is negative i.e. of a minus, the 
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function gives a value of -1. It is essential that the 
value is 2 times larger than the conventional minmod 
limiter function. By using the AF ( I # J) , it is only 
required to define the slope function as: 
5 G(I,J)=AF(I,J)/Ax ...(25) 

An expression Sp(I,J)/Ax means a forward 
difference from the corresponding one of the 
definitions of the above equations (23), while an 
expression Sm(I,J)/Ax means a backward difference from 

10 the corresponding one of the definitions of the same. 
If the terms Sp(I,J) and Sm(I,J) have different signs, 
the sum (sgn(SpI,J) + sgn(Sm(I,J)) becomes equal to 
zero, and hence the whole limiter function AF(I,J) and 
therefore the slope function G(I,J) becomes equal to 

15 zero. If the terms do not have different signs, the 

slope function G(I,J) becomes equal to a value twice as 
large as the smaller one of the absolute values of the 
forward difference and the backward difference, since 
the sum (sgn(SpI,J) + sgn(Sm(I,J)) becomes equal to 2 

20 or -2. 

. As described above, according to the interpolation 
method of the present embodiment, with respect to a 
cell of interest on the one-dimensional structured grid, 
the slope of the function H(I,J) is set to zero if the 
25 signs of the forward difference and the backward 

difference of the mapping F on the cell are different 
from each other, and to a value twice as large as the 
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smaller one of the absolute values of the differences 
if the signs are the same, the value F(X,T) at the 
center of the cell of interest is determined as the 
linear function, and based on the linear function thus 
5 determined, values of the function H(I,J) on a partial 
region of the one-dimensional real region determined by 
the cell are defined. 

FIG. 1 is a flowchart showing a process of 
advection calculation using the above described 

10 interpolation method according to the first embodiment. 

First, in a step SI, initialization is performed, 
e.g. by generating a grid, determining time step 
intervals, and then initializing the initial shape of a 
function F ( • , 1) . 

15 In the following step S2, the computation of the 

above -described equation Q(-,J) (see the equations (20) 
and (21)) is performed. At this time, as shown in the 
above equations (20) and (21), the interpolation 
function H according to the present embodiment is used. 

20 Then, in a step S3, based on data obtained by the 
equation Q(*,J), the computation of the function 
F(-,J+1) is carried out (see the above equation (22)). 

Thereafter, time is updated in a step S4 . More 
specifically, the problem being solved now is a time 

25 evolution equation expressed by the above equation (18) , 
and according to the conventional practice of numerical 
solution of an ordinary time evolution equation, the 
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value F(-,J) of the function F at time T=AtJ is 
determined from the value F(-,J-1), and then time is 
changed from T=At*(J-l) to T=At*J. 

Then, in a step S5, it is determined whether or 
5 not the present step is a final step, and if the 

present step is not a final step, the process returns 
to the step S2 , whereby the present computation is 
performed until it is determined that the present step 
is a final step. 

10 The use of the interpolation method according to 

the present embodiment provides two separate regions, 
one defined by a function with a slope of zero, and the 
other defined by a function with a sharp slope, which 
makes it possible to obtain very good results as a 

15 function describing the shape of a two-phase 

incompressible fluid. FIG. 2 shows a result of the 
computation. The present computation is carried out 
under the same conditions corresponding to those of the 
conventional computations results of which are shown in 

20 FIGS. 7 to 9, and the shape corresponding to the 

initial value and that after 1000 steps are depicted in 
FIG. 2. It will be understood that the result of 
advection calculation carried out by the interpolation 
method according to the present embodiment shown in FIG. 

25 2 maintains a very sharp rectangular shape compared 
with those shown in FIGS. 7 to 9 . 

Next, a description will be given of an apparatus 
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for carrying out the interpolation method according to 
the present embodiment. 

FIG. 3 is a perspective view of a computer system 
for executing advection calculation using the 
5 interpolation method shown in FIG. 1, and FIG. 4 is a 
block diagram showing -the arrangement of essential 
components of the computer system. 

As shown in FIGS. 3 and 4, the computer system for 
carrying out the interpolation method is comprised of a 

10 processor including a plurality of CPU's 116 and a 

memory 118, a storage device including a hard disk 117 
and a floppy (registered trademark) disk 110, an input 
section including a keyboard 115 and a mouse 114, and 
an output display section including a display 112. The 

15 processor and other components are accommodated in a 
main body 111 . 

This computer system loads a program designed 
according to the flowchart shown in FIG. 1 and encoded, 
into the memory 118, reserves storage for necessary 

20 computation, performs predetermined arithmetic 

operations based on a shape of a region, coordinates of 
grid points, and physical quantities, which have been 
input by an appropriate method, and writes values of 
the physical quantities obtained on the region by the 

25 arithmetic operations into the hard disk 117 or the 

like for storage and/or causes the same to be displayed 
on the display 112 . 
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The use of the interpolation method according to 
the present embodiment makes it possible to form a 
function on the real axis from values given on the grid 
points though the values include discontinuities. 
5 Therefore, even when the present invention is applied 
to a usage different from the present embodiment, e.g. 
data extraction in signal processing, or one- 
dimensional image processing, it is possible to express 
a sharp shape. 

10 Next, a description will be given of an 

interpolation method according to a second embodiment 
of the present invention. In the second embodiment, 
the interpolation algorithm for a one-dimensional grid, 
described hereinabove as to the first embodiment, is 

15 extended to one for a two-dimensional grid. 

The equation dealt with in the present embodiment 
is a partial differential equation in a two-dimensional 
space, defined by: 

5F(x,y,t)/3 t = <?Cx(x,y, t)F(x,y, t) / <?x 

20 + <9Cy(x,y,t)F(x,y,t)/<?y ...(26) 

The problem of converting this function to a 
difference equation for computation will now be 
considered. Suppose a function F(X,Y,T) as a 
discretized form of F(x,y,t). The arguments X, Y, in 

25 the function F cooperatively represent a discrete 
spacial grid point and the argument T in the same 
represents a discrete time. 



First, it is assumed that an interpolation 
function H[I, J,T] (x,y) according to the present 
embodiment has already been obtained as the following 
piecewise linear function: 
5 H[I,J,T] (x,y)=F(I, J,T)+Gx(I, J # T)x+Gy(I, J,T)y .,.(27) 

With this assumption, a description will be given 
of a method of numerically solving the above equation 
(26) by finding a solution thereto by finite difference 
approximation. A method of determining coefficients Gx 

10 (I,J,T) and Gy(I/J,T) of the interpolation function 
H[I, J/T] (x,y) of the present embodiment will be 
described hereinafter. As described later, the method 
of determining the coefficients Gx(I,J,T) and Gy(I,J / T) 
characterizes the interpolation algorithm according to 

15 the present embodiment. It is essential to the present 
embodiment that the coefficients Gx(I,J,T) and 
Gy(I,J,T) are determined by using neighboring 
information of geometric properties of the grid. 
In short, when there are provided a two- 

20 dimensional structured grid formed on a two-dimensional 
real region, and a function F(X,Y,T) on the two- 
dimensional structured grid, the function being defined 
through definition of a value thereof at the center of 
each cell within the two-dimensional structured grid, 

25 the function F(X,Y,T) is defined as an interpolation 
function H[I, J,T] (x,y) on the two-dimensional real 
region^ by using the function F(X,Y,T) and the geometry 
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of the grid. 

If the function Cx(I+l/2, J,T) is positive, there 

is given the following definition: 

Qx (1+1/2, J,T) :=Cx( 1+1/2, J,T)H[I, J,T] 
( Ax/2-Cx(I, J,T) At/2, -Cy (I, J,T) At/2) ...(28) 
If the function Cx ( 1+1/2 , J, T) is not positive, 

there is given the following definition: 

Qx( 1+1/2, J,T) :=Cx( 1+1/2, J,T)H[I+1, J,T] 
(- Ax/2+Cx(I, J,T) At/2, -Cy(I, J,T) At/2) ...(29) 
If the function Cy ( I , J+l/2 , T) is positive, there 

is given the following definition: 

Qy(I, J+l/2, J,T) :=Cy(I, J+1/2,T)H[I, J,T] 
(-Cx(I, J,T) At/2, Ay/2-Cy(I, J,T) At/2) ...(30) 
If the function Cy ( I , J+l/2 , T) is not positive, 

there is given the following definition: 

Qy(I, J+l/2, J,T) :=Cy(I / J+l/2,T)H[I / J+l / T] 
(-Cx(I,J,T) At/2,-Ay/2-Cy(I,J,T) At/2) ...(31) 
From the above definitions, there is obtained the 

following equation: 

F(I, J,T+1) =F(I, J,T) + At (Qx( 1+1/2, J, T) 
-Qx (1-1/2, J,T) ) / Ax+ At (Qy(I, J+l/2, T) 
-Qy(I,J-l/2,T))/Ay ...(32) 

On the other hand, a slope function G(I,J,T) 
representative of the slope of the piecewise linear 
function is defined in the following manner: 

First, there are introduced the following 
definitions : 
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Sxp(I,J,T) :=F(I+1, J,T) -F(I, J,T) 

Sxm(I,J,T) :=F(I, J,T)-F(I-1, J,T) 

Syp(I,J,T) :=F(I, J+1,T)-F(I, J,T) 

Sym(I,J,T) :=F(I, J,T) -F(I, J-1,T) ...(33) 

Further, there are introduced limiter functions 

one dimension for the present embodiment, defined by 
Sx(I, J)=(sgn(Sxp(I, J,T) ) +sgn (Sxm (I , J, T) ) ) 
Min( |Sxp(I, J,T) | , | Sxm (I, J, T) | ) 
Sy(I, J)=(sgn(Syp(I, J,T) ) +sgn (Sym ( I , J, T) ) ) 
Min( |Syp(I, J,T) | , |Sym(I, J,T) | ) ...(34) 
Then, there are provided the following 

definitions : 

Vppmin=Min(F(I, J+1,T) -F(I, J,T) , F ( 1+1 , J+l , T) 
-F(I, J,T) , - £ ,F(I + 1, J,T) -F(I, J,T) ) 
Vmpmin=Min(F(I-l, J+1,T) -F(I, J,T) ,F(I, J+1,T) 
-F(I,J,T) ,F(I-1, J,T)-F(I, J,T) ,- 6 ) 
Vpmmin=Min(- £ ,F (1+1, J,T) -F(I, J,T) , F(I, J-1,T) 
-F(I, J,T) ,F(I+1, J-1,T) -F(I, J,T) ) 

Vmmmin=Min(F (1-1, J,T)-F(I, J,T) ,- £ , F (1-1 , J-l , T) 
-F(I, J,T) ,F(I, J-1,T) -F(I, J,T) ) ...(35) 
Vppmax=Max(F(I, J+1,T) -F(I, J,T) , F ( 1+1, J+l , T) 
-F(I, J,T) , - £ ,F(I + 1, J,T) -F(I, J,T) ) 
Vmpmax=Max(F (1-1, J+1,T) -F(I, J,T) , F(I, J+1,T) 
-F(I, J,T) ,F(I-1, J,T) -F(I, J,T) , - £ ) 
Vpmmax=Max(- £ ,F(I+1, J,T) -F(I, J,T) ,F(I, J-1,T) 
-F(I, J,T) ,F(I+1, J-1,T)-F(I, J,T) ) 
Vmmmax=Max(F(I-l, J,T) -F(I, J,T) ,- £ , F (1-1 , J-l, T) 
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-F(I, J,T) ,F(I, J-l^J-Ftl, J,T) ) ..,(36) 

There are also provided the following definitions: 

hpp(I, J,T)=Sx(I, J,T) /2+Sx(I, J,T) 12 

hmp(I, J.Tl^Sxtl, J,T) /2+Sx(I / J,T) 12 
5 hpm(I, J,T)=Sx(I, J,T) /2-Sx(I, J,T) 12 

hiran(I, J,T)=-Sx(I, J,T) /2-Sx(I / J,T) 12 ...(37) 

Further, there is performed the following 
calculation: 

Assuming, for the time being, that an equation 
10 phipp=phimp=phipm=phiram=l holds, if hpp>Vppmax holds, 

an equation of phipp=Vppmax/hpp is newly defined, while 
if hpp<Vppmin holds, an equation of phipp=Vppmin/hpp is 
newly defined; 

Similarly, if hpm>Vpmmax holds, an equation of 
15 phipm=Vpmmax/hpm is newly defined, while if hpm<Vpmmin 
holds, an equation of phipm=Vpmmin/hpm is newly 
defined; 

Also, if hmp>Vmpmax holds, an equation of 
phimp=Vmpmax/hmp is newly defined, while if hmp<Vmpmin 

20 holds, an equation of phimp=Vmpmin/hmp is newly 
defined; and 

If hmm>Vmmmax holds, an equation of 
phiiruna=Vnunmax/hmm is newly defined, while if hmm<Vmmmin 
holds, an equation of phiinm=Vmmmin/hmm is newly defined. 

25 By using the thus obtained values phipp, phimp, 

phipm, phimm, the above-mentioned coefficients 
Gx(I,J,T) and Gy(I,J,T) are determined by the following 
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Gx(I, J,T) 

=Min(phipp,phipm,phimp,phimm) Sx(I, J,T) /Ax 
Gy(I, J,T) 

5 =Min { phipp , phipm , phimp , phimm )Sy(I,J,T) / Ay ...(38) 

Thereafter, the interpolation function 
H(I, J,T) (x,y) within the grid (I, J) is defined as: 

H(I,J,T) (x,y)=F(I f J,T)+Gx(I,J / T)x+Gy(I / J,T)y ...(39) 
The interpolation method according to the present 
10 embodiment can be summarized as follows: 

A cell of interest is set to a cell A, the cell A 
having a first side extending in an x direction, a 
second side extending in the x direction and being 
opposite to the first side, a third side extending in a 
15 y direction, and a fourth side extending in the y 

direction and being opposite to the third side, and 
values twice as large as values of one-sided 
differences of the above described function F between 
the center of the cell A and respective centers of four 
20 cells adjacent to the cell A on the first side, the 
second side, the third side, and the fourth side are 
defined as a first-sided difference (DFxmin) , a second- 
sided difference (DFxmax) , a third-sided difference 
(DFymin) , and a fourth-sided difference (DFymax) , 
25 respectively. If the first-sided difference and the 
second-sided difference have different signs, an x- 
direction difference is set to zero, whereas if the 
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first-sided difference and the second-sided difference 
have the same sign, the x-direction difference is set 
to the smaller one of the absolute values of the first- 
sided difference and the second-sided difference. If 
5 the third-sided difference and the fourth-sided 
difference have different signs, a y-direction 
difference is set to zero, while if the third-sided 
difference and the fourth-sided difference have the 
same sign, the y-direction difference is set to the 

10 smaller one of the absolute values of the third-sided 
difference and the fourth-sided difference. There is 
formed an interpolation function candidate which has 
slopes defined by the x-direction difference and the y- 
direction difference, respectively, and has a value F0 

15 of the above-mentioned function F at the center of the 
cell A, and this candidate is set to a plane B. If a 
value of the plane B at a first grid point at a 
location of intersection of the first side and the 
third side of the cell A is larger than the value of 

20 the center of the cell A, the plane B is modified by 
multiplying the x-direction difference and the y- 
direction difference by a largest constant not more 
than 1 such that the value of the plane B at the first 
grid point does not exceed any of the values of the 

25 function F at respective centers of three cells having 
the first grid point in common except for the cell A. 
If the value of the plane B at the first grid point is 
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smaller than the value of the center of the cell A, the 
plane B is modified by multiplying the x-direction 
difference and the y-direction difference by a largest 
constant not more than 1 such that the value of the 
5 plane B at the first grid point does not fall below any 
of the values of the function F at the respective 
centers of the three cells having the grid point in 
common except for the cell A. On the plane B thus 
obtained, the same operation as carried out as to the 

10 first grid point is carried out as to a second grid 

point at a location of intersection of the first side 
and the fourth side of the cell A, a third grid point 
at a location of intersection of the second side and 
the third side, and a fourth grid point at a location 

15 of intersection of the second side and the fourth side, 
to thereby change the slope of the plane B, and the 
resulting plane is defined as the interpolation 
function H within the cell A. This algorithm is 
equivalent to the algorithm described hereinbefore. 

2 0 FIG. 5 is a flowchart showing a process of 

advection calculation using the above described 
interpolation method according to the second embodiment. 

First, in a step Sll, initialization is performed, 
e.g. by generating a grid, determining time step 

25 intervals, and then initializing the initial shape of a 
function F ( • , • , 1) . 

In the following step S12, the computation of the 
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above-described equation Q(-,-,J) (see the equations 
(28) to (31) ) is performed. At this time, as shown in 
the above equations (28) and (29), the interpolation 
function H according to the present embodiment is used. 
5 Then, in a step S13, based on data obtained by the 
equation Q(',',J), the computation of the function 
F(-,',J+1) is carried out (see the above equation (32)). 

Thereafter, time is updated in a step S14 . More 
specifically, the problem being solved now is a time 

10 evolution equation expressed by the above equation (26) , 
and according to the conventional practice of numerical 
solution of an ordinary time evolution equation, the 
value F(-,-,J) of the function F at time T=AtJ is 
determined from the value F(',-,J-1), and then time is 

15 changed from T=At* (J-l) to T=At*J. 

Then, in a step S15, it is determined whether or 
not the present step is a final step, and if the 
present step is not a final step, the process returns 
to the step S12, whereby the present computation is 

20 performed until it is determined that the present step 
is a final step. 

The use of the interpolation method according to 
the present embodiment provides two separate regions, 
one defined by a function with a slope of zero, and the 

25 other defined by a function with a sharp slope, which 
makes it possible to obtain very good results as a 
function describing the shape of a two-phase 
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incompressible fluid . FIGS. 6A and 6B show results of 
the computation. The present computation is carried 
out under the same conditions corresponding to those of 
the conventional computations results of which are 
5 shown in FIGS. 10A and 10B, and the shape corresponding 
to the initial value and that after 1000 steps are 
depicted in FIGS. 10A and 10B, respectively. It will 
be understood that the result of advection calculation 
carried out by the interpolation method according to 
10 the present embodiment shown in FIG. 6B maintains a 

very sharp rectangular shape compared with those shown 
in FIG. 10B. 

The computer for executing the advection 
calculation using the interpolation function according 
15 to the present embodiment is identical to those shown 
in FIGS. 3 and 4. 

Next, a description will be given of an 
interpolation method according to a third embodiment of 
the present invention. In the third embodiment, the 
20 interpolation algorithm for a two-dimensional grid, 
described hereinabove as to the second embodiment is 
extended to one for a three-dimensional grid. 

The equation dealt with in the present embodiment 
is a partial differential equation in a three- 
25 dimensional space, defined by: 

dF(x,y,z,t)/<9 t=c?Cx(x,y,z,t)F(x,y,z,t)/(9x 
+ dCy (x,y, z, t)F(x,y, z, t) / dy 



+ c?Cz(x,y,z,t)F(x,y,z,t)/dz ... (40) 

The problem of converting this function to a 
difference equation for computation will now be 
considered. Suppose a function F(X,Y,Z,T) as a 
5 discretized form of F(x,y,z,t). The arguments X, Y, Z 
in the function F cooperatively represent a discrete 
spacial grid point and the argument T in the same 
represents a discrete time. 

First, it is assumed that an interpolation 
10 function H[I, J, K, T] (x,y, z) according to the present 

embodiment has already been obtained as the following 
piecewise linear function: 

H[I, J,K,T] (x,y / 2)=F(I,J,K / T)+Gx(I / J f K,T)x 

+Gy(I, J,K,T)y+Gz (I, J,K,T) z ...(41) 
15 With this assumption, a description will be given 

of a method of numerically solving the above equation 
(40) by finding a solution thereto by finite difference 
approximation. A method of determining coefficients Gx 
(I,J,K,T), Gy(I,J,K,T) , and Gz (I, J, K, T) of the 
20 interpolation function H [ I , J, K, T] (x, y, z ) of the present 
embodiment will be described hereinafter. As described 
later, the method of determining the coefficients 
Gx(I,J,K,T), Gy(I,J,K,T), and Gz(I,J,K,T) characterizes 
the interpolation algorithm according to the present 
25 embodiment. It is essential to the present embodiment 
that the coefficients Gx(I,J,K,T), Gy ( I , J, K, T) , and 
Gz(I,J,K,T) are determined by using neighboring 
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information of geometric properties of the grid. 

In short, when there are provided a three- 
dimensional structured grid formed on a three- 
dimensional real region, and a function F(X,Y,Z,T) on 
the three-dimensional structured grid, the function 
being defined through definition of a value thereof at 
the center of each cell within the three-dimensional 
structured grid at a given time T, the function 
F(X,Y,Z,T) is defined as an interpolation function 
H[I,J,K, T](x,y,z) on the three-dimensional real region 
by extending and interpolating the function F(X,Y,Z,T) 
through utilization of the function F(X,Y,Z,T) and the 
geometry of the grid. 

If the function Cx (1+1/2, J, K,T) is positive, 
there is given the following definition: 

Qx( 1+1/2, J,K,T) :=Cx (1+1/2 , J,K,T)H[I, J,K,T] 

( Ax/2-Cx(I, J,K,T) At/2, -Cy(I, J,K,T) At/2, 

-Cz(I, J,K,T) At/2) ...(42) 

If the function Cx (1+1/2 , J, K, T) is not positive, 
there is given the following definition: 

Qx( 1+1/2, J,K,T) :=Cx( 1+1/2, J, K, T) H [ 1+1 , J,K,T] 
(- Ax/2+Cx(I, J,K,T) At/2,-Cy(I, J,K,T) At/2, 
-Cz (I, J,K,T) At/2) ...(43) 

If the function Cy ( I , J+l/2 , K, T) is positive, there 
is given the following definition: 

Qy(I, J+l/2, K,T) :=Cy(I, J+l/2 , K, T) H [I, J,K,T] 
(-Cx(I, J,K,T) At/2, Ay/2-Cy (I, J,K,T) At/2, 



-Cz(I,J,K,T) At/2) ...(44) 

If the function Cy (I, J+l/2 ,K,T) is not positive, 
there is given the following definition: 

Qy (I, J+1/2,K,T) :=Cy(I, J+l/2 , K, T) H [I , J+1,K,T] 
(-Cx(I, J,K,T) At/2 / -Ay/2-Cy(I / J / K / T) At/2, 
-Cz(I, J,K,T) At/2) ...(45) 

Further, if the function Cz (I, J, K+l/2 , T) is not 
positive, there is given the following definition: 
Qz (I, J, K+l/2, T) :=Cz (I , J, K+l/2 , T) H [I , J, K, T] 
(-Cx(I, J,K,T) At/2, -Cy (I, J,K,T) At/2, 
A z/2-Cz (I, J, K, T) At/2) ...(46) 

If the function Cy ( I , J, K+l/2 , T) is not positive, 
there is given the following definition: 

Qz (I, J, K+l/2, T) :=Cz (I, J, K+l/2 , T) H [ I , J,K+1,T] 
(-Cx(I, J,K,T) At/2, -Cy (I, J,K,T) At/2, 
-Az/2-Cz (I, J,K,T) At/2) ...(47) 

From the above definitions, there is obtained the 
following equation : 

F(I, J,K,T+1)=F(I, J,K,T)+At (Qx( 1 + 1/2, J,K,T) 
-Qx( 1-1/2, J,K,T) ) / Ax+At (Qy(I, J+l/2, K,T) 
-Qy(I, J-1/2,K,T) ) / Ay+At (Qz (I, J, K+l/2 , T) 
-Qz(I, J,K-1/2,T) ) / Az ...(48) 

On the other hand, a slope function G(I,J,K,T) 
representative of the slope of the piecewise linear 
function is defined in the following manner: 

First, there are introduced the following 
definitions : 
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Sxp(I, J,K,T) :=F(I+1, J,K,T) -F(I, J,K,T) 

Sxm(I, J,K,T) :=F(I, J,K,T) -F (I-l , J, K, T) 

Syp(I, J,K,T) :=F(I,J+1,K,T)-F(I,J,K,T) 

Sym(I, J,K,T) : =F ( I , J, K, T) -F ( I , J-l , K, T) 

Szp(I, J,K,T) :=F(I, J,K+1,T) -F(I, J,K,T) 

Szm(I, J,K,T) :=F(I,J,K,T)-F(I,J,K-1,T) ...(49) 

Further, there are introduced limiter functions in 

one dimension for the present embodiment, defined by: 
Sx(I, J,K)=(sgn(Sxp(I, J,K,T) ) +sgn (Sxm (I , J, K, T) ) ) 
Min( |Sxp(I, J,K,T) | , |Sxm(I, J,K,T) | ) 
Sy(I, J,K)=(sgn(Syp(I, J,K,T) )+sgn(Sym(I, J,K,T) ) ) 
Min( |Syp(I, J,K,T) | , |Sym(I, J,K,T) | ) 
Sz(I, J,K)=(sgn(Szp(I, J,K,T) ) +sgn (Szm ( I , J, K, T) ) ) 
Min( |Szp(I, J,K,T) | , | Szm ( I , J, K, T) | ) ...(50) 
Then, there are provided the following 

definitions : 

Vpppmin=Min(F(I, J+1,K,T) -F(I, J,K,T) , F ( 1+1 , J+l , K, T) 
-F(I, J,K,T) , - £ ,F(I+1, J,K,T)-F(I, J,K,T) , 
F(I, J+1,K+1,T) -F(I, J,K,T) ,F(I+1, J+1,K+1,T) 
-F(I, J,K,T) ,F(I, J,K+1,T)-F(I, J,K,T) , 
F(I+1, J,K+1,T) -F(I, J,K,T) ) 

Vmppmin=Min(F(I-l, J+1,K,T) - (I, J,K,T) , F (I , J+l , K, T) 
-F(I, J,K,T) ,F(I-1, J,K,T) -F(I, J,K,T) , - £ , 
F(I-1, J+1,K+1,T) -F(I, J,K,T) ,F(I, J+l, K+1,T) 
-F(I, J,K,T) ,F(I-1,J,K+1,T)-F(I, J,K,T) , 
F(I, J,K+1,T) -F(I, J,K,T) ) 

Vpmpmin=Min ( - £ , F ( 1+1 , J, K, T) -F ( I, J, K, T) , 



48 

F(I,J-1,K,T)-F(I,J,K,T) ,F(I+1, J-1,K,T) -F(I,J,K,T) 
F(I, J,K+1,T) -F(I, J,K,T) ,F(I+1, J,K+1,T) -F(I, J,K,T) 
F(I,J-1,K+1,T)-F(I,J,K,T) ,F(I+1, J-1,K+1,T) 
-F(I, J,K,T) ) 

Vmmpmin=Min(F(I-l, J,K,T) -F ( I , J, K, T) , - £ , 

F(I-1, J-1,K,T)-F(I, J,K,T) ,F(I, J-1,K,T) -F(I, J,K,T) 

F(I-1, J,K+1,T) -F(I, J,K,T) ,F(I, J,K+1,T) -F ( I , J,K,T) 

F(I-1, J-1,K+1,T) -F(I, J,K,T) , F ( I , J-l , K+l , T) 

-F(I, J,K,T) ) 

Vppnunin=Min(F(I, J+1,K,T) - (I, J,K,T) , F ( 1 + 1 , J+l , K, T) 
-F(I, J,K,T) ,- £ ,F(I + 1, J,K,T)-F(I, J,K,T) , 
F(I, J+1,K-1,T) -F(I, J,K,T) ,F(I+1, J+1,K-1,T) 
-F(I, J,K,T) ,F(I,J,K-1,T)-F(I,J,K,T) , 
F(I+1, J,K-1,T)-F(I, J,K,T) ) 

Vmpmmin=Min(F(I-l, J+1,K,T) - (I, J,K f T) , F ( I , J+l , K, T) 
-F(I, J,K,T) ,F(I-1, J,K,T) -F(I, J,K,T) , - £ , 
F(I-1,J+1,K-1,T) -F(I, J,K,T) ,F(I, J+1,K-1,T) 
-F(I, J,K,T) ,F(I-1, J,K-1,T)-F(I, J,K,T) , 
F(I, J,K-1,T)-F(I, J,K,T) ) 

Vpmmmin=Min ( - £ , F ( 1+1 , J, K, T) -F(I, J,K,T) , 

F(I, J-1,K,T)-F(I, J,K,T) ,F(I+1, J-1,K,T) -F(I, J,K,T) 

F(I, J,K-1,T)-F(I, J,K,T) ,F(I+1, J,K-1,T) -F(I, J,K,T) 

F(I, J-1,K-1 # T) -F(I, J,K,T) ,F(I+1, J-1,K-1,T) 

-F(I, J,K,T) ) 

Virararanin=Min(F(I-l / J,K,T) -F (I,J,K,T),-£, 

F(I-1, J-1,K,T)-F(I, J,K,T) ,F(I, J-1,K,T) -F(I, J,K,T) 

F ( 1-1 , J, K-l , T) -F ( I , J, K, T) ,F(I, J,K-1,T) -F(I,J,K,T) 
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P(I-1,J-1,K-1,T)-F(I,J,K,T) ,F(I, J-1,K-1,T) 
-F(I, J,K,T) ) ...(51) 

Vpppmax=Max (F ( I , J+l , K, T) - ( I , J, K, T) , F ( 1+1 , J+l , K, T) 
-F(I, J,K,T) , - £ ,F(I+1, J,K,T)-F(I, J,K,T) , 
F(I, J+l, K+l, T) -F(I, J,K,T) ,F(I+1, J+1,K+1,T) 
-F(I, J,K,T) ,F(I, J,K+1,T)-F(I, J,K,T) , 
F(I+1, J,K+1,T)-F(I, J,K,T) ) 

Vmppmax=Max(F(I-l, J+l f K,T) - (I, J,K,T) , F (I, J+l , K, T) 
-F(I,J,K,T),F(M,J,K,T)-F(I,J,K,T),-£, 
F(I-1, J+1,K+1,T) -F(I, J,K,T) , F(I, J+1,K+1,T) 
-F(I,J,K,T) ,F(I-1, J,K+1,T)-F(I, J,K,T) , 
F(I, J,K+1,T) -F(I, J,K,T) ) 

Vpmpmax=Max ( - £ , F ( 1+1 , J, K, T) -F(I, J,K,T) , 

F(I, J-1,K,T)-F(I, J,K,T) ,F(I+1, J-l , K, T) -F ( I , J,K,T) , 

F ( I , J, K+l , T) -F(I, J,K,T) ,F(I+1, J,K+1,T) -F(I, J,K,T) , 

F (I, J-1,K+1,T) -F(I, J,K,T) , F(I + 1, J-l, K+l ,T) 

-F(I, J,K,T) ) 

Vmmpmax=Max(F(I-l, J,K,T) -F(I, J,K,T) , - £ , 

F(I-1, J-1,K,T)-F(I, J,K,T) ,F(I, J-1,K,T) -F(I,J,K,T), 

F(I-1, J,K+1,T)-F(I, J,K,T) ,F(I, J,K+1,T) -F(I, J,K,T) , 

F(I-1, J-l, K+l, T) -F(I, J,K,T) , F(I, J-l, K+l, T) 

-F(I, J,K,T) ) 

Vppitimax=Max(F(I, J+1,K,T)-(I, J,K,T) , F ( 1 + 1 , J+l , K, T) 
-F(I, J,K,T) ,- £ ,F(I + 1, J,K,T) -F(I, J,K,T) , 
F(I, J+1,K-1,T) -F(I, J,K,T) ,F(I+1, J+1,K-1,T) 
-F(I, J,K,T) ,F(I, J,K-1,T)-F(I, J,K,T) , 
F(I+1, J,K-1,T)-F(I, J,K,T) ) 
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Vmpirumax=Max(F(I-l, J+1,K,T) - (I, J,K,T) , F ( I , J+l, K, T) 
-F(I, J.K.T) ,F(I-1, J,K,T)-F(I, J,K,T) ,- 6 , 
F(I-1, J+1,K-1,T) -F(I, J,K,T) ,F(I, J+1,K-1,T) 
-F(I, J,K,T) ,F(I-1, J,K-1,T)-F(I, J,K,T) , 
F(I, J,K-1,T) -F(I, J,K,T) ) 

Vpmmmax=Max ( - £ , F (1+1 , J, K, T) -F(I, J,K,T) , 

F(I, J-1,K,T) -F(I, J,K,T) ,F(I+1, J-1,K,T) -F(I, J,K,T) , 

F(I,J,K-1,T) -F(I, J,K,T) ,F(I+1, J, K-l , T) -F ( I , J,K,T) , 

F(I, J-1,K-1,T) -F(I, J,K,T) ,F(I+1, J-1,K-1,T) 

-F(I, J,K,T) ) 

Vmmmmax=Max(F(I-l, J,K,T) -F (I, J,K,T) , - £ , 
F(I-1, J-1,K,T) -F(I, J,K,T) ,F(I, J-l, K, T) -F (I, J,K,T) , 
F(I-1, J,K-1,T) -F(I, J,K,T) ,F(I, J, K-l , T) -F ( I , J,K,T) , 
F(I-1, J-1,K-1,T)-F(I, J,K,T) ,F(I, J-l, K-l, T) 
-F(I,J,K,T)) ...(52) 

There are also provided the following definitions: 
hppp(I, J,K,T) 

=Sx(I, J,K,T) /2+Sx(I, J,K,T) /2+Sz ( I , J, K, T) /2 
hmpp(I, J,K,T) 

=-Sx(I, J,K,T) /2+Sx(I, J,K,T) /2+Sz(I, J,K,T) /2 
hpmp(I, J,K,T) 

=Sx(I, J,K,T) /2-Sx(I, J,K,T) /2+Sz(I, J,K,T) /2 
hmmp (I, J, K, T) 

=-Sx(I, J,K,T) /2-Sx(I, J,K,T) /2+Sz ( I , J, K, T) 12 
hppm(I, J,K,T) 

=Sx(I, J,K,T) /2+Sxd, J,K,T) /2-Sz (I , J, K, T) /2 
hmpmfl, J,K,T) 
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= -Sx(I, J,K,T) /2+Sx(I, J,K,T) /2-Sz ( I , J, K, T) /2 
hpmmd, J,K,T) 

=Sx(I, J,K,T) /2-Sx(I, J,K,T) /2-Sz ( I , J, K, T) /2 
hmmm ( I , J, K, T) 

5 =-Sx(I, J,K,T) /2-Sx(I,J / K,T) /2-Sz (I, J,K,T) /2 ...(53) 

Further, there is performed the following 
calculation: 

Assuming, for the time being, that an equation 
phippp=phimpp=phipmp=phiimp=phippm=phimpm=phipmm=phiiranm 
10 =1 holds, 

if hppp>Vpppmax holds, an equation of 
phippp=Vpppmax/hppp is newly defined, while if 
hppp<Vpppmin holds, an equation of phippp=Vpppmin/hppp 
is newly defined; 
15 Similarly, if hpmp>Vpmpmax holds, an equation of 

phi pmp = Vpmpmax / hpmp is newly defined, while if 
hpmp<Vpmpmin holds, an equation of phipmp=Vpmmin/hpmp 
is newly defined; 

Also, if hmpp>Vmppmax holds, an equation of 
2 0 phimpp=Vmppmax/hmpp is newly defined, while if 

hmpp<Vmppmin holds, an equation of phimpp=Vmppmin/hmpp 
is newly defined; 

If hmmp>Vmmpmax holds, an equation of 
phimmp=Vmmpmax/hmmp is newly defined, while if 
25 hmmp<Vmmpmin holds, an equation of phimmp=Vmmpmin/hnuup 
is newly defined; 

If hppm>Vppmmax holds, an equation of 
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phippm=Vppmmax/hppm is newly defined, while if 
hppm<Vppmmin holds, an equation of phippm=Vppmmin/hppm 
is newly defined; 

If hpmm>Vpmmmax holds, an equation of 
phipmm=Vpmmmax/hpmm is newly defined, while if 
hpmm<Vpmmmmin holds, an equation of phipim=Vpmmmin/hpinm 
is newly defined; 

If hmpm>Vmpmmax holds, an equation of 
ph i mpm= Vmpmmax / hmpm is newly defined, while if 
hmpm< Vmpmm i n holds, an equation of ph i mpm= Vmpmm in/ hmpm 
is newly defined; and 

If hirmun>Vmmmmax holds, an equation of 
phimim=Vmmmmax/hmmp is newly defined, while if 
hmmm<Vmmmmin holds, an equation of phinmm\=Vmmmmin/hmmm 
is newly defined. 

By using the thus obtained values phippp, phimpp, 
phipmp, phimmp, phippm, phimpm, phipmm, phimmm, the 
above -mentioned coefficients Gx(I,J,K,T), Gy(I,J,K,T), 
and Gz(I,J,K,T) are determined by the following 
equations : 

Gx ( I , J, K , T) =Min (phippp , phimpp , phipmp , phimmp , phippm, 
phimpm, phipmm, phimmm) Sx ( I , J, K, T) / Ax 

Gy ( I , J , K , T) =Min (phippp , phimpp , phipmp , phimmp , phippm , 
phimpm, phipmm, phimmm) Sy (I, J, K, T) / Ay 

Gz (I , J, K, T) =Min (phippp , phimpp , phipmp , phimmp , phippm 
phimpm, phipmm, phimmm) Sz ( I , J, K, T) / A z 

Thereafter, the interpolation function 
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H(I, J,K,T) (x,y, z) within the grid (I,J,K) is defined 
as : 

H[I, J,K,T] (x,y,z)=F(I,J / K,T)+Gx(I / J / K,T)x 
+Gy (I, J,K,T)y+Gz (I, J,K,T) z ...(54) 

The interpolation method according to the present 
embodiment can be summarized as follows: 

A cell of interest is set to a cell A, the cell A 
having a first side (face, in this three-dimensional 
grid) extending in an x direction, a second side 
extending in the x direction and being opposite to the 
first side, a third side extending in a y direction, a 
fourth side extending in the y direction and being 
opposite to the third side, a fifth side in a z 
direction, and a sixth side extending in the z 
direction and being opposite to the fifth side, and 
values twice as large as values of one-sided 
differences of the above described function F between 
the center of the cell A and respective centers of six 
cells adjacent to the cell A on the first side, the 
second side, the third side, the fourth side, the fifth 
side, and the sixth side are defined as a first-sided 
difference (DFxmin) , a second-sided difference (DFxmax) , 
a third-sided difference (DFymin) , a fourth-sided 
difference (DFymax) , a fifth-sided difference (DFzmax) , 
and a sixth-sided difference (DFzmin) , respectively. 
If the first-sided difference and the second-sided 
difference have different signs, an x-direction 



54 

difference is set to zero, whereas if the first-sided 
difference and the second-sided difference have the 
same sign, the x-direction difference is set to the 
smaller one of the absolute values of the first-sided 
5 difference and the second-sided difference. If the 

third-sided difference and the fourth-sided difference 
have different signs, a y-direction difference is set 
to zero, while if the third-sided difference and the 
fourth-sided difference have the same sign, the y- 

10 direction difference is set to the smaller one of the 
absolute values of the third-sided difference and the 
fourth-sided difference. If the fifth-sided difference 
and the sixth-sided difference have different signs, a 
z-direction difference is set to zero, while if the 

15 fifth-sided difference and the sixth-sided difference 
have the same sign, the z-direction difference is set 
to the smaller one of absolute values of the fifth- 
sided difference and the sixth-sided difference. There 
is formed an interpolation function candidate which has 

2 0 slopes defined by the x-direction difference, the y- 
direction difference, and the z-direction difference, 
respectively, and has a value F0 of the above-mentioned 
function F at the center of the cell A, and this 
candidate is set to a plane B. If a value of the plane 

25 B at a first grid point at a location of intersection 
of the first side, the third side, and the fifth side 
of the cell A is larger than the value of the center of 
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the cell A, the plane B is modified by multiplying the 
x-direction difference, the y-direction difference, and 
the z-direction difference by a largest constant not 
more than 1 such that the value of the plane B at the 
5 first grid point does not exceed any of the values of 
the function F at respective centers of seven cells 
having the first grid point in common except for the 
cell A. If the value of the plane B at the first grid 
point is smaller than the value of the center of the 

10 cell A, the plane B is modified by multiplying the x- 
direction difference, the y-direction difference, and 
the z-direction difference by a largest constant not 
more than 1 such that the value of the plane B at the 
first grid point does not fall below any of the values 

15 of the function f at the respective centers of the 

seven cells having the grid point in common except for 
the cell A. 

On the plane B as the candidate interpolation 
function thus obtained, the same operation as carried 

2 0 out as to the first grid point is carried out as to a 
second grid point at a location of intersection of the 
first side, the third side, and the sixth side, a third 
grid point at a location of intersection of the first 
side, the fourth side, and the fifth side, a fourth 

25 grid point at a location of intersection of the first 

side, the fourth side, and the sixth side, a fifth grid 
point at a location of intersection of the second side, 
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the third side, and the fifth side, and a sixth grid 
point at a location of intersection of the second side, 
the third side, and the sixth side, a seventh grid 
point at a location of intersection of the second side, 
5 the fourth side, and the fifth side, and an eighth grid 
point at a location of intersection of the second side, 
the fourth side, and the sixth side, to thereby change 
the slope of the plane B, and the resulting plane is 
defined as the interpolation function H within the cell 

10 A. This algorithm is equivalent to the algorithm 
described hereinbefore. 

A flowchart of an advection calculation using the 
interpolation method of the present embodiment is the 
same as shown in FIG. 5. Further, the interpolation 

15 method according to the third embodiment can be 

realized by the same computer described hereinbefore 
with reference to FIGS. 3 and 4. In the present 
embodiment as well, it is possible to carry out the 
advection calculation while maintaining a sharp shape 

20 of the function as is the case with the interpolation 
method for a two-dimensional space according to the 
above described second embodiment. 

If the interpolation methods according to the 
second and third embodiments are used, it is possible 

25 to form a function in a two-dimensional or three- 
dimensional space based on values given to grid points 
of the space although the values include 
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discontinuities. By virtue of this feature, even when 
the interpolation method according to the present 
invention is applied to other uses than the above 
described embodiments, e.g. reproduction of image data, 
5 detailing data of a digital two-dimensional or three- 
dimensional object, it is possible to express a sharp 
shape in the same manner. 

Furthermore, the present invention is not limited 
to the apparatuses according to the above described 
10 embodiments but may be either applied to a system 

composed of a plurality of apparatuses or to a single 
apparatus . 

It is to be understood that the object of the 
present invention may be accomplished by supplying a 

15 system or an apparatus with a storage medium in which a 
program code of software which realizes the functions 
of any of the above described embodiments, and causing 
a computer (CPU or MPU) of the system or apparatus to 
read out and execute the program code stored in the 

2 0 storage medium. 

In this case, the program code itself read from 
the storage medium realizes the functions of any of the 
above described embodiments, and hence the storage 
medium on which the program code is stored constitutes 

25 the present invention. 

Examples of the storage medium for supplying the 
program code include a floppy (registered trademark) 
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disk, a hard disk, an optical disk, a magneto optical 
disk, a CD-ROM, a CD-R, a CD-RW, a DVD-ROM, a DVD-RAM, 
a DVD-RW, a DVD+RW, a magnetic tape, a nonvolatile 
memory card, and a ROM. Downloading via a network can 
also be utilized. 

Further, it is to be understood that the functions 
of any of the above described embodiments may be 
accomplished not only by executing a program code read 
out by a computer, but also by causing an OS (operating 
system) or the like which operates on the computer to 
perform a part or all of the actual operations based on 
instructions of the program code. 

Further, it is to be understood that the functions 
of any of the above described embodiments may be 
accomplished by writing a program code read out from 
the storage medium into a memory provided on an 
expansion board inserted into a computer or in an 
expansion unit connected to the computer and then 
causing a CPU or the like provided in the expansion 
board or the expansion unit to perform a part or all of 
the actual operations based on instructions of the 
program code . 



